#---
#FUNCTIONS: SIMULATE NONHOMOGENOUS POISSON PROCESS FOR VISITING TIMES
#---

#nhpp: generates a sequence of event times using a Gamma distribution intensity function
#note: uses inversion method algorithm
#m - number of events
#T - maximum time
#s - shape parameter
#r - rate parameter
nhpp <- function(m, T, s, r) 
{
	Nevents <- rpois(1, m)
	if(Nevents == 0)
	{
		result <- NULL
	}

	else
	{
		result <- rgamma(Nevents, shape = s, rate = r)
		while(any(result > T))
		{
			result[result > T] <- rgamma(sum(result > T), shape = s, rate = r)
		}
		result <- sort(result)
	}
	
	return(result)
}
	
#sim.nhpp.data: simulates 'n' sequences of event times based on 'nhpp' function
#n - number of simulations
#m - number of events
#tuning - parameter to adjust probabilities of missing a visit (between 0 and 1)
#T - maximum time
#s - shape parameter
#r - rate parameter
sim.nhpp.data <- function(n, m, tuning = 0.25, T, s, r)
{
	#simulate visit times based on arrival probability
	tt <- 1:T
	u <- matrix(runif(length(tt) * n), nrow = n)
	visit <- matrix(rep(tt, each = n), nrow = n)
	arr.prob <- 1.05 - (tt/T)^tuning
	visit.time <- visit * t(apply(u, 1, function(u) ifelse(u <= arr.prob, 1, NA)))

	#simulate event times for each individual
	sim.data <- NULL
	sim.data.exact <- NULL
	for (i in 1:n)
	{
		event.time <- nhpp(m, T, s, r)
		if(!is.null(event.time))
		{
			sim.data.exact <- rbind(sim.data.exact, data.frame(ID = i, TIME = event.time))
			visiting <- sort(unique(visit.time[i,visit.time[i,] != 0 & visit.time[i,] != T]))
			ind.count <- table(1 + findInterval(event.time, visiting))
			temp.data <- data.frame(ID = i, L = c(0, visiting), R = c(visiting, T), M = 0)
			temp.data[as.numeric(names(ind.count)), 4] <- ind.count
			sim.data <- rbind(sim.data, temp.data) 
		}
	}

	#censored event times
	temp.data.int <- data.frame(ID = unique(sim.data$ID), 
		CT = sapply(split(sim.data$R, sim.data$ID), max))
	sim.data.int <- merge(sim.data, temp.data.int)
 
	#exact event times
	sim.data.exact <- merge(sim.data.exact, temp.data.int)
	temp.data.exact <- data.frame(ID = unique(sim.data.exact$ID), 
		N = sapply(split(sim.data.exact$TIME, sim.data.exact$ID), length))
	sim.data.exact <- merge(sim.data.exact, temp.data.exact)

	#return censored and exact times
	return(list(int = sim.data.int, exact = sim.data.exact))
}
